Vortex matter in layered superconductors without Josephson coupling: 
numerical simulations within a mean-field approach 
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We study vortex matter in layered superconductors in the limit of zero Josephson coupling. 
The long range of the interaction between pancake vortices in the c-direction allows us to employ 
a mean-field method: all attractive inter-layer interactions are reduced to an effective substrate 
potential, which pancakes experience in addition to the same-layer pancake repulsion. We perform 
numerical simulations of this mean-field model using two independent numerical implementations 
with different simulation methods (Monte-Carlo sampling and Langevin molecular dynamics). The 
substrate potential is updated self-consistently from the averaged pancake density. Depending on 
temperature, this potential converges to a periodic profile (crystal) or vanishes (liquid). We compute 
thermodynamic properties of the system, such as the melting line, the instability line of the crystal, 
and the entropy jump across the melting transition. The simulation results are in good agreement 
with approximate analytical calculations. 

PACS numbers: 74.60.Ge, 74.80.Dm 



I. INTRODUCTION 

The vortex state in type-II superconductors is 
a complex physical system. Within the layered 
high-temperature materials, such as BiaSraCaCuaOs 
(BSCCO) or YBa 2 Cu 3 7 (YBCO), vortex lines ea 
understood as wiggling stacks of pancake vortices .Be 
The thermodynamic properties of the vortex state are 
determined by the interaction between pancake vortices. 
There are two mechanisms of pancake interaction: (i) 
electromagnetic interaction and (ii) Josephson coupling. 
The electromagnetic interaction is mediated by super- 
currents circulating around each pancake, whereas the 
Josephson coupling results from the energy cost due to a 
phase shift between the superconducting order parame- 
ters in the neighboring layers. 

To understand the phase diagram of high-temperature 
superconductors and in particular the melting line of the 
vortex lattice,ErQ we need to gain an insight into the be- 
havior of vortex matter under a variety of experimental 
conditions. In moderately anisotropic materials, such as 
YBCO, the short-range Josephson coupling is the domi- 
nant inter-layer interaction- [ a^ r the- | vortices are well de- 
scribed as elastic strings J3'DE3EJc3li3 In very anisotropic 
materials on the other hand, such as BSCCO, the Joseph- 
son coupling is weak, and the long-range electromagnetic 
interaction between the pancakes should be taken into 
account. In this paper we consider very anisotropic ma- 
terials in the absence of Josephson coupling, and neglect 
pinning. Even after keeping only the electromagnetic 



coupling, the problem remains a challenging one, due to 
the long range of the interactions: the energy of elec- 
tromagnetic interaction between two pancakes depends 
logarithmically on the separation along the layers and 
decays exponentially with the number of layers between 
the pancakes. More specifically, the interaction is repul- 
sive between pancakes in the same layer, and attractive 
between pancakes in different layers, and the decay length 
of the exponential dependence is the London penetration 
depth A, which is typically 100 times larger than the layer 
spacing s. Approximately, this system has beep investi- 
gated within the density-functional theoryO't£lE3 

For a numerical investigation of the system, one can 
in principle simulate directly a stack of two-dimensional 
(2D) pancake systems taking into account all of the inter- 
layer interactions. However, the computational challenge 
is that the interlayer attraction between pancakes ex- 
tends over a range of 2X/s ~ 100— 150 layers. In addition, 
realistic simulation of the melting transition requires at 
least several hundred point vortices per layer. So far, di- 
rect numerical investigations have been performed only 
on small systems using aJjavit. 10 layers and of the or- 



der of 100 vorticesOHHEErEilH This is not sufficient to 
describe realistically the vortex state in BSCCO. With 
today's computational resources, it is not feasible to per- 
form realistic direct three-dimensional (3D) simulations 
of this system because the necessary computational effort 
grows quadratically with the number of layers. 

Fortunately, one can benefit from the long range of the 
interlayer coupling. As the interlayer force on a pancake 
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FIG. 2: Left: Snapshot of the pancakes (visualized by 
spheres) placed onto the substrate potential (visualized as 
a surface) at B = B\ and t = 1/59 ~ 0.017 just below the 
melting transition to demonstrate how the substrate potential 
constrains the pancake motion. Right: The phase diagram we 
have computed using the substrate method. 



FIG. 1: Schematic representation of the substrate model. 
(A) The pancake positions p n (x) in each layer n, are (B) av- 
eraged over the layers in order to obtain the averaged pancake 
density p(x). From the average pancake density we compute 
(G) the substrate potential VmfOe), which is smeared over a 
length of the order of A. The vortex lattice spacing is an and 
s is the layer spacing. 



is the result of a sum of a large number (~ 2X/s) of 
small contributions, it can be calculated by a mean- field 
approach. The exact value of this force is determined by 
the instantaneous pancake densities in the large number 
of layers. In the crystal state the instantaneous den- 
sity can be decomposed into the average density, which 
is a periodic function of the in-plane coordinates, and a 
fluctuating contribution. In the mean-field approach to 
the interlayer interactions, one replaces the instantaneous 
densities in the other layers by the average density. This 
approach gives a quantitatively correct description of the 
system, because due to the law of large numbers, the ne- 
glected force from the fluctuating densities is typically 
smaller than the average interlayer force by the factor 
~ y/s/X <C 1. The calculation then takes the form of in- 
dependent layers, with the pancakes, in each layer feeling 
an effective "substrate potential" r3 This substrate po- 
tential is the cumulative affect of the attraction of pan- 
cakes in all other layers as illustrated in Fig. Pan- 
cakes within one layer interact directly with each other, 
whereas the interaction with pancakes in other layers is 
mediated via the substrate potential. Thus, each layer is 
treated individually, until a new substrate potential can 
be computed. This process is iterated, until the substrate 
has converged to a steady solution. In this paper, we 
present the first numerical implementations of this sub- 
strate model and show results which we compare with 
the semi-analytic approximations given in Ref. p3[ . 

We summarize this work in Fig. |^. On the left the 
central idea is visualized: pancakes experience attractive 
inter-layer interactions through the substrate potential 



which stabilizes the pancake crystal. On the right we 
show the computed melting line separating a 3D pancake 
vortex lattice from decoupled 2D liquids. We express 
magnetic induction in units of B\ = <E>n/A 2 , where $o 
is the magnetic flux quantum, such that the pancake- 
spacing in a triangular lattice is an = \J 2/V3 A rs 1.07A 
at B = B\. We use a dimensionless temperature, t, which 
is the ratio of the thermal energy k^T to the prefactor 
2sen of the logarithmic pancake-pancake interaction, 
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where en = Qq/^it figX 2 ), /in is the vacuum permeability, 
and s is the layer spacing. This allows to compare our 
results with outcomi^iMMrL 2D one-component Coulomb 
plasma simulations ,r'Hr 5 H 26 l where frequently T = 1/t is 
used to express temperatures. At low fields, the electro- 
magnetic attraction of range A S> s between pancakes 
in different layers stabilizes the 3D pancake-vortex lat- 
tice. Increasing the magnetic field decreases the rela- 
tive strength of the inter-layer coupling. At high fields, 
B 3> B\, the long-range repulsive interaction within the 
layer dominates, and the 3D pancake lattice melts at a 
temperature close to the 2D melting temperature. 

In Sec. H we describe the substrate model in detail, 
including three different methods for t he effi cient compu- 
tation of the substrate potential (Sec. II C ). The results, 
including the equilibrium phase dia gram , are shown in 
Sec. Ill before we conclude in Sec. IV. The appendix 
gives a derivation of the correlation correction to the free 
energy, and shows that our mean-field approach should 
be accurate to order s/X. 
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II. MEAN-FIELD APPROACH (SUBSTRATE 
MODEL) 

A. The mean-field inter-layer coupling 

The in-layer energy £ 111 and the inter-layer energy £ lnter 
of a system of electromagnetically interacting pancakes 
in a layered superconductor is, respectively, 

£ in = £ = £ \ U(Rj - R™/ , 0) (2) 



and 



f inter = 1 £ Yu(R]-Kf,,n-n'). (3) 

Indices n and n' count over layers, and j and j' over pan- 
cakes in the layers, R" is the (2D) position of pancake 
j in layer n, and U(R, n) is the coupling energy for two 
pancakes separated by a vector (R, z), where z = ns, 
with s being the layer spacing. The z-axis is chosen per- 
pendicular to the layers. n 
The in-layer pancake interaction.! is 



(4) 



JL f°° w exp(-rVA) 
2A ./„ r' 



and the inter-layer interaction (n ^ 0) is 
1 ' ns\ . f L 



U(r,n) — — |<xj) 



A 



( ns \ 



In | - 

r 



(5) 



dr 



, exp(- v /r' 2 + (ns) 2 /A) 



Using 



we rewrite 



(6) 



n — n'). 



£ i„ter = i 53 / d 2 rd 2 r' p„(r)p„,(r')^(r - r', 

We separate pancake density fluctuations from the layer- 
average density 



p(r) = (p„(r)), 
Pn(r) = p(r) +<fy>„(r), 



(8) 
(9) 



and obtain from (f7j) 

f totet = 1 ^ / d 2 r dV[/(r - r', n - n') (10) 
x [p{v)p{v') + 2p{v')5 Pn {Y) + 5p n (r)5 Pn ,(r')] . 



Because the difference n — n' in the last sum extends 
over a very large number of layers (~ A/s), a typical 
value of the sum J2 n ' U(r — r',n — n')2p(r') is larger 
than a typical value of the sum J2 n ' U(r — r' 7 n~ n')8p n i 
by the factor ~ yj A/ s. Again, the law of large numbers 
allows us to neglect the last term in (10) leading to the 
mean-field description of the interlayer interactions. A 
more precise justification is given in the Appendix, where 
the free energy correction due to the correlation term is 
shown to be smaller than the mean-field free energy by 
the factor s/A. 

Separating the pancake density into the average value 
P = (Pn(r)) and a modulating part, we can split the to- 
tal magnetic coupling energy into two parts, each with a 
quite different meaning. The part containing the average 
density does not depend on temperature and formally di- 
verges due to the logarithmic term in (||) ■ This divergence 
exactly compensates a similar divergence in the in-plane 
energy. Within the mean-field approach the part of the 
coupling energy sensitive to density variations is finite 
only in the crystal state. In the liquid state it vanishes. 

For the mean-field interlayer energy £ MF we obtain 
from @ 

£ MF = Y,E™ F (11) 

n 

= i /d 2 rdVf/(r-r',n-n') 

x[p(r)p(v') + 2p(r')6p n (r)} 
= \ E [ d2r Wr)p(r) +J2 [ d2r VMF(r)6p n (r) 

(12) 

- - ^E / d * r ^MF(r)p(r) + £ fd 2 r V MF (r)p n (r). 

(13) 

The last term describes fluctuations in the fixed substrate 
potential Vmf, 



dV 



with 



= Jd 2 r'U(r - r')p(r') 
= (U*p)(r) 



U{v) = Y,U{v,n). 

njtO 



p(r') (14) 

(15) 
(16) 

(17) 



U(r) is the interaction potential of a pancake separated 
by r from a stack of pancakes minus the interaction q£ 
the (missing) pancake in the same layer and is given byu 



«(r) = 2e oS K (j) -U{r,0) 



(18) 
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with Kq (x) being a modified Bessel function of the second 
kind. Ignoring terms of the order of s/2A, the pancake- 
pancake repulsion (||) simplifies to 



U(r,0) = 2e sln 



L 



(19) 



In our calculations we find it useful to use the form in 
Fourier space ,EZI 



Uicv : 4 "' !l> ^ A- 2 + g 2 ~7 

A- 2 



= —Alters 



<7 2 (A- 2 + g 2 ) 



(20) 
(21) 



B. Algorithm 

In principle, the substrate model can be implemented 
as follows: 

1. Assume initial pancake densities /9 n (r), for example 
a hexagonal lattice in each layer n. 

2. Average the pancake density p n (r) over all layers 
to obtain p(r), (||). 

3. Compute the substrate potential Vmf (r) , (|l6|), by 
convoluting the substrate interaction kernel U{r), 
(|l8|), with the average pancake density p(r) 



V MF (r) = (U*p)(r). 



(22) 



4. For each layer n compute the pancake distribu- 
tion p n (r) using Monte Carlo or Langevin dynamics 
simulations. The total energy for layer n contains 
the direct pancake-pancake interaction within the 
layer® 



(23) 



and the relevant interaction with pancakes in other 
layers via the substrate potential 



E 



MF 



d 2 r V MF (r)p(r) + d 2 r% F (r)p n (r) 



_Ecoup 



-£ coup + $>mf(R"), 



(24) 



E coup is constant for a given p(r) and can there- 
fore be ignored within the Monte Carlo/Langevin 
simulation as it only shifts the energy scale. 

5. Go to 0, until Vmf ( or p) has converged. 



Since the substrate potential Vmf in step [| is the same 
for all layers, we can compute p n (r) for many Langevin- 
dynamics time-steps (or Monte-Carlo sweeps) rather 
than many layers. Therefore, in order to obtain the av- 
eraged pancake density p(r) in step 0^ we average over 
time-steps (or sweeps) computed in one layer rather than 
averaging over layers. 

Using the substrate potential, we reduce the solution of 
the 3D problem to performing one 2D simulation in the 
presence of the iteratively refined substrate potential. 



C. Numerical implementation 

We exploit the convolution theorpin and compute the 
substrate potential in Fourier-space 



Vmf(i-) ^ (U*p)(r) 



(25) 



/i2 
W(q)p(q)e X p(iq-r) (26) 

using the analytical Fourier transform U(q) as given in 
(pl|), and the numerically computed 

p(q) = y"d 2 r p(r)exp(-iq-r). (27) 

This has two advantages: firstly, we do not cut off the 
interaction kernel U within the simulation cell as would 
be the case in the real-space convolution. Secondly, this 
is numerically more efficient than performing the convo- 
lution ( |l5| ) directly. 

We have used three different methods for computing 
Vmf(f) numerically. 



1. The full method 

The "full method" computes the substrate potential 
Vmf using the full spectrum p(q) of Fourier-components 
of the average pancake density p(r) as shown in (|2q). 
In our simulations we use a resolution of « 100 2 grid- 
cells per pancake in order to compute p(r) as an av- 
erage over time-steps/sweeps. This results in recipro- 
cal lattice vectors up to magnitudes of « lOOQo, where 
Qo = 47r/(V3a ), because |Q max |/Qo ~ 2vr/(AxQo) ~ 
aO/Ax « 100. The necessary discrete Fourier-transform 
of p(r), and the inverse transform of Vmf (q) = W(q)/o(q) 
can be done efficiently rtising an implementation of the 
Fast Fourier TransformEa 

We pre-compute the substrate potential VmfO") on a 
mesh and interpolate subsequently for intermediate pan- 
cake positions while performing Langevin dynamics in 
the fixed substrate. We compute a new substrate every 
200,000 time-steps. It is important to average over so 
many time-steps to reduce density fluctuations (due to 
poor statistics) in the pancake histogram, which would 
result in a deformed substrate potential. 
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Note that p(r) and p(q) are discretized out of numeri- 
cal necessity to compute a histogram but not for concep- 
tual reasons. 



2. The Fourier-filtered method 

The average density p(r) should be a periodic function, 
which can be represented by a discrete set of Fourier 
components. Therefore, the second method uses only a 
subset Q M of the Fourier components q to represent p(q) 



,FF 



(q) = (27r) 2 5> QM( 5 2 (q-Q, 



(28) 



which we determine from the maxima of the structure 
factor and 



1 



PQ 



L X Ly 

" 3 



exp(-iR^-Q) 



(29) 



with L x L y being the area of the simulation cell. We 
average over a set of configurations c of pancake positions 
(either sweeps or time-steps) to compute pq. 

Using p FF (r) = (2tt)~ 2 f d 2 q p FF (q) exp(iq-r) to 
present p(r), we Fourier-filter p(r), and keep only the 
relevant components for the computation of the periodic 
substrate. We can write 



d 2 q 



W(qV FF (q)exp( iq .r) (30) 



(2n) 2 

^W(Q A1 )pq M exp(iQ M -r) 



(0) „^ PQ» exp(zQ M -r) 

- -^osj: (32) 



This is equivalent to using the full-method, but setting 
p(Q) - if Q g {CU. 

The advantage of the Fourier-filtered method is that 
we need to average over less iterations before we can 
compute a new pancake density, and subsequently a new 
substrate, because the substrate is per construction peri- 
odic. Using the Fourier-filtered method we use 500 time- 
steps/sweeps for each substrate iteration. 

It turns out that it is not necessary to take the average 
( p9| ) over different configurations but it is sufficient to use 
just one configuration (i.e. one time step or sweep): 



PQ 



1 



L X L, 



■ exp(— iRj-Q). 



(33) 



Nevertheless, we run a simulation for 500 time- 
steps/sweeps with the same fixed substrate potential to 
reduce re-computation of pq, and to give the pancakes 
some time to explore the system with a new substrate 
potential. 
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FIG. 3: Set of Q M -vectors up to third order (i.e. three 
"shells" around the origin) in reciprocal space used in the 
small-harmonics Fourier-filtered method to compute pQ . 
Due to the reality of p(v) we have p(q) = p(— q) and it is 
therefore sufficient to compute only half of the 36 coefficients 



3. The small-harmonics (Fourier- filtered) approximation 

In addition to Fourier-filtering p(r) we can speed up 
the computation further because close to the melting 
temperature, pq M decays quickly for higher-order Q M due 
to the Debye- Waller factor. We can estimate the reduc- 
tion of pQ due to the Debye- Waller factor 



exp 



(u 2 )Q 2 



exp 



1 (u 2 ) 16tt 2 Q< 



4 a 2 



(34) 



(31) where Qq = 47r/(v / 3ao) and (u 2 ) is the mean-squared 



fluctuation displacement. Depending on (u ) we can 
ignore all jo(Q M ) with |Q M | > Q'. For all but the 
smallest fields, we find close to t he me lting transition 
(u 2 )/al « 0.02 - 0.03 (see Section |HI E| ), and it is suffi- 
cient to include up to 3 rd -order vectors Q M in the sum- 
mation in ( |32] ) as shown in Fig. [| 

For the small-harmonics Fourier-filtered method it is 
more efficient to evaluate (|32| ) for each pancake posi- 
tion occurring in the Langevin/Monte-Carlo simulation 
rather than pre-computing Vmf on a mesh. 

We demonstrate the equivalence of the full and the 
Fourier-filtered method for the determination of the in- 



stability line in Sec. [II A, and we compare with the 
small-harmonics Fourier-filtered method in Sec. Ill C. 



D. Monte Carlo and Langevin dynamics 
simulations 



We have two independent implementations of the 
small-harmonics Fourier-filtered method: AEK has writ- 
ten a Monte Carlo simulation that is based on energy 
evaluations, and HF has implemented a Langevin dynam- 
ics simulation based on force calculations. The results of 
both implementations agree perfectly. 
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FIG. 4: Convergence to pancake lattice at B — B\ and 
t — 1/59 » 0.017. Top: pancake histogram n(x,yo) taken 
along y = yo- Bottom: substrate potential Vmf0e,2/o)- 
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FIG. 5: Convergence to pancake liquid. As Fig. ^ but at 
t — 1/50 = 0.02 above the melting temperature at B = B\. 



We follom. standard vortex-state simulation 
techniques J^ju including periodic boundary condi- 
tions for the in-plane interactions. We use a .smooth 
cut-off for the vortex in-plane interactions.ESEil For 
the Langevin dynamics simulations, we compute the 
substrate f orces numerically from the pre-computed 
mesh (Sec. II CI) for the full method and the Fourier- 



filtered method. For the small-harmonics Fourier-filtered 
method we use the analytical derivative of (j||). The 
Monte Carlo simulations were only implemented with 
the small-harmonics method. If not stated otherwise we 
use a system with 1020 pancakes. 



III. RESULTS 



A. Time convergence of the substrate potential 



As described in Sec. [IB, we start each run with a 



hexagonal pancake distribution corresponding to zero 
temperature. Fig. [I| shows results for the Fourier-filtered 
method at B = B\ and at a temperature t = 1/59 « 
0.017. The top plot shows a one-dimensional slice of 
the 2D pancake histogram n(x,yo) taken along x at 
y = yo. The histogram relates to the pancake density via 
n(x,y) = p(x,y)AxAy where Ax and Ay are the spac- 
ings of the grid used to create the histogram. For the th 
substrate-iteration we set the histogram to have narrow 
and high peaks at the pancake equilibrium positions cor- 
responding to delta-peaks in a zero-temperature pancake 
density p(r). Based on this initial pancake distribution, 
we compute the substrate potential, Vmf(i")j f° r tne nrst 
substrate-iteration, of which a one-dimensional slice at 
y = yo is shown in the lower part of Fig. |4| Using this 
substrate potential, we run the Langevin dynamics sim- 
ulation for 200,000 steps which results in the histogram 
for iteration 1 as shown in the upper plot of Fig. 0. 



Based on these data, we compute the substrate potential 
for iteration 2. We iterate the substrate re-computation 
until the substrate potential has reached a steady state 
(after typically than 10 substrate-iterations). The fig- 
ure demonstrates that the system converges quickly to a 
pancake solid at this temperature below melting. 

The dotted line in the lower part of Fig. [| shows a 
comparison substrate potential for iteration 10 computed 
using the full method. While the amplitude and width 
of the wells (and thus the resulting force) are virtually 
identical to the Fourier-filtered data, the magnitude of 
the substrate from the full method varies slightly. This 
is due to (long wavelength) density fluctuations in the 
histogram data and reduces further if one uses more time- 
steps for each substrate iteration. 

Fig. H shows data for B = B\ and a higher tempera- 
ture t = 1/50 = 0.02 which is above the melting temper- 
ature. Here, the pancake distribution broadens and con- 
sequently the substrate potential flattens quickly within 
the first few substrate iterations. Eventually, the system 
has become a disordered liquid with an approximately 
constant pancake density and the substrate is virtually 
flat, as shown for iteration 10. We conclude that for this 
temperature and magnetic field the pancake lattice is un- 
stable to melting into a pancake liquid. 

For Fig. H and || we have used 200,000 time-steps for 
each substrate-iteration in order to be able to compare 
the full and the Fourier-filtered method, but it would 
be sufficient to use much less time-steps per substrate- 
iteration for the Fourier-filtered methods. For pro- 
duction purposes, we use the small-harmonics Fourier- 
filtered met hod an d update the substrate every 500 time- 
steps (Sec. [I C 2 ). Although more substrate- iterations 
than with the full method are required before the sys- 
tem reaches a steady state, the small-harmonics Fourier- 
filtered approach is more efficient. The full method and 
both Fourier-filtered methods find that at B = B\ the 
pancake lattice becomes unstable for 0.017 < t < 0.018. 



7 



1000 2000 



0.0300 


I i 
- i 


i 1 i 


w 0.0275 


_ V . • 


* * : 


O 

cu 






™ 0.0250 






H 

CD 






0.0225 






0.0200 


: i 


i , F 



FIG. 6: Finite-size investigation of instability temperature 
at B = 0AB x . 
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FIG. 7: Example of a hysteresis loop obtained by heating a 
crystal (circles) and cooling a liquid (squares) for B — 0AB\. 
Each point was equilibrated for 9 • 10 Monte-Carlo steps. 
The crystal melts at t = 0.027, while the liquid freezes at 
t = 0.0247. Rhombs represent results obtained by simula- 
tions starting from the same intermediate defective configu- 
ration with ridcf ~ 0.2, this configuration melts at t = 0.0262, 
which we take as an estimate for the thermodynamic melting 
temperature. The insets show dependencies of the defect con- 
centration on the Monte Carlo step at the temperatures where 
the intermediate configuration melts and the liquid configu- 
ration freezes. See also [B2f. 



Fig. |6| shows how the instability temperature varies as 
a function of system size. For small numbers of vortices, 
iV" p , the temperature oscillates slightly and for larger sys- 
tems it becomes constant. Most importantly, there is 
no general trend visible although the data ranges from 



N p = 90 to ./Vp = 1512. This insensitivity to the system 
size demonstrates the local nature of the melting transi- 
tion at this field. 



B. Hysteresis loop 

Rather than starting from a hexagonal crystal for every 
temperature, a better approach to determine the insta- 
bility temperature is to subsequently increase t until the 
system melts. We also find that by starting from a liq- 
uid configuration and lowering t, the system jumps into 
the crystal state at a certain freezing temperature. We 
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3 rd order Q 
* 20 lh order Q 

□ infinite lattice summation (3 rd order Q) 




0.005 



0.015 



0.02 
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FIG. 8: Top: The substrate curvature a s at B = B\. Shown 
is our numerical solution of the Fourier-filtered method us- 
ing Q- vectors up to 20 th order (stars), and using up to 3 rd 
order (circles). We have also shown results for Q ma x = 3Qo 
using an infinite lattice summation for the in-layer interac- 
tion (squares). Bottom: The pancake fluctuation width (u 2 ) 
(stars). For comparise«, we also show the results of a sim- 
ple SCHA calculations (solid line) of the softening of the 
substrate potential (not including thermal softening effects of 
the 2D lattice). The results are close to our numerical data 
at low temperatures, but as the melting point is approached 
there are extra (anharmonic) fluctuations in the simulation 
data for (u 2 ), resulting in an even softer substrate potential 



expect the true melting temperature to lie within the in- 
stability and the freezing temperature. Such a hysteretic 
run is shown in Fig. 0. 



In order to estimate the thermodynamic melting tem- 
perature at which the free energy of the solid and the liq- 
uid phase cross, we proceed as follows. Firstly, we store 
a vortex configuration taken from 2D melting-transition 
simulations. We chose a configuration from a time- 
step/sweep where the system was previously a solid but 
just starts melting, i.e. the defect density starts shooting 
up and the structure factor peaks start decaying. This 
vortex configuration is "intermediate" between a solid 
and a liquid. Secondly, we start the computation from 
this intermediate configuration (IC) for every tempera- 
ture. The results for the IC simulations are shown in 
Fig. ^ (rhombs). We use the temperature at which the 
IC melts as a best approximation to the melting temper- 
ature of the physical 3D pancake- vortex lattice. 
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C. Temperature dependence of substrate curvature 
and pancake fluctuation width 

We can quantify the strength of the substrate potential 
with 



(35) 



This is the curvature of the potential evaluated at de- 
viations Rj — R° from the equilibrium lattice positions 
and averaged over pancake positions R 7 . The second 
derivative can be taken analytically from (p2L). 

Fig. H shows in the upper plot how a s varies with tem- 
perature. The solid line is an analytical prediction from 
treating the substrate softening due to thermal fluctua- 
tions within the self-consistent harmonic approximation 
(SCHA) £3 All other data are simulation results from the 
Fourier-filtered method. The stars show a s computed 
using the small-harmonics Fourier-filtered method with 
Q M -vectors up to 20 th order. For low temperatures the 
data nearly coincide with the SCHA-solution. Close to 
melting the SCHA-a s is larger than the numerical result. 
Therefore the simulations give a softer substrate and the 
lattice has larger thermal displacements. This difference 
could be due to the inadequacies of the SCHA which 
does not include the thermal softening of the 2D lattice. 
The more complex two-vertex self-consistent harmonic, 
approximation (2VSCHA) does include these effects.^ 
The circles in Fig. |8| show results using Q M -vectors up to 
3 rd order, as shown in Fig. ||. Close to the transition 
from solid to liquid around t w 0.0175 these data agree 
perfectly with the higher-order data. At lower temper- 
atures the 3 rd order results deviate from the 20 th order 
because (u 2 ) becomes smaller in the Debye- Waller factor 
(j34|). However, as long as we are interested in temper- 
atures close to the transition, the 3 rd order approach is 
sufficient. 

The square boxes are computed using the 3 rd order 
approach, but instead of smoothly reducing the pan- 
cake interactionEil at a distance of « 7an, we use an 
infinite lattice summation technique for the logarithmic 
interaction.^! This demonstrates that it is sufficient to 
use a (smooth) cut-off for the in-layer pancake interac- 
tions. 

We compute the average pancake fluctuation width 
(it 2 ) by fitting to a distribution where each pancake is 
normally smeared around its equilibrium position R° 



27TO- 2 ^ 

j 



exp 



R?| 2 



2cr 2 



(36) 



The Fourier transform of p(r) is 



p(q) 



exp 



i;2 ^) 2 no exp(^-)^ ( 5 2 (q-Q M ). 
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FIG. 9: Phase diagram of the electro-magnetically coupled 
3D pancake system. Numerically computed instability line 
(black circles on dashed line) in comparison with the insta- 
bility line from the 2VSCHA (dash-dotted line). Also shown 
is an semi-analytical estimate for the melting line from Ref. 
|23| , to be compared with our numerical estimate tic (crosses). 
We have shaded the solid phase underneath the melting in 
gray. The melting temperature t„ = 0.007 of a 2D system is 
shown by a dotted line. 



The Fourier components p(Q M ) have the Debye- Waller 
factor as an envelope, and by fitting a Gaussian to it, we 
can determine (it 2 ) = 2er 2 . 

The lower part of Fig. || shows computed values for 



We express 



in units of ai and it increases 



from at zero temperature towards 0.028 close to the 
transition, which corresponds to a Lindemann number 
of ps 0.168 at B = B\. In agreement with an over- 
estimation of a s by the SCHA, (it 2 ) is underestimated in 
comparison with the numerical results close to the melt- 
ing transition. 



D. Phase diagram 



As demonstrated in Sec. Ill A, we can determine for 



each parameter pair (B,T) whether the pancake system 
remains a 3D pancake lattice, or whether it is unstable 
towards the liquid phase which consists of decoupled 2D 
liquids. (This is sometimes called a pancake gas, even 
though there are still very strong in-plane correlations 
in the decoupled layers. In the absence of Josephson 
coupling, a line-like liquid regime is expected only at ex- 
tremely small magnetic fields.E-3) 

We probe par amete r space in the i?-T-plane as de- 
scribed in Sec. |IIIA| and compute an instability line 
for the phase diagram of the system, which is shown in 
Fig. H (circles on dashed line). We also show an esti- 
mate of the instability line that has been computed using 
the two-vertex self-consistent harmonic approximate 
(2VSCHA, dash-dotted) for the substrate model 
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j v per pancake across the melting 
The inset shows the same data 



FIG. 10: Top: Latent heat L p 
transition as a function of field, 
on a reduced scale. Bottom: Jump in inter-layer coupling 
energy, AI7 coup = E coup , normalized by latent heat. 




0.0 -; 



FIG. 11: Entropy jump AS P per pancake as a function of 
field B. Inset: Pancake fluctuation width (it 2 ) normalized to 



at the melting point (the Lindemann number squared). 



Since in this work we explicitly compute the pancake 
positions without using approximations (within the sub- 
strate model), we expect our result to be more accurate 
than the 2VSCHA. It can be seen that the 2VSCHA 
slightly overestimates the temperature for the instabil- 
ity line. 

Our numerical estimates for the melting points at cer- 
tain fields (see Section 1IIB) are shown as the crosses 
in Fig. H Also shown as the solid line is the melting 
line calculated semi-analytically in Ref. ^3|. In this work 
the melting temperature was estimated by comparing ap- 
proximate free energies F = U — TS for the solid and 
liquid phases. The solid free energy was calculated from 
the SCHA, which gives a variational upper bound on the 
free energy. The liquid free energy was taken from earlier 
simulations of a single layer J£j i.e. it was assumed that 
the layers are completely uncoupled in the liquid state. 
Remarkably, our melting points from simulations lie on 
top of the semi-analytic line (to within our error bars). 

For increasing fields B, the substrate becomes weaker 
and weaker and the melting temperature drops. In the 
limit of B — ► oo we recover a 2D system with logarithmic 
interactions for which melting has been estimateoOi£3'E£l 
to occur at T^f w 140 ± 10 44> t 2 ^ « 0.007, which is 
consistent with our results. 

At low fields the pancake stacks are widely separated 
and interact only weakly with each other. In this limit 
the system melts below the,-ejvaporation transition of an 
isolated stack of pancakeaffia at T = 4 o t = 0.25. In 
agreement with this, we find that the instability line ap- 
proaches t w 0.25 for B — > (see Fig. 0). 



E. Latent heat and jump in entropy 

We compute the latent heat per pancake, L p , by taking 
the difference of internal energy between the solid and the 



liquid phase at the melting temperature T r 
1 



L n — tt~ (^liquid — [/solid) (37) 

(i^d + £ coup ))- (38) 



liquid 



The internal energy U of one layer in the solid phase 
consists of the in-plane energy E lu (23) and the inter- 
layer coupling energy E coup , whereas iJ cou P = in the 
liquid phase in our model. In order to compute E coup 
for the solid phase, we use ( p"2| ) where the second sum 
vanishes due to the definition of 8p n : 



£cou P= 1 |d 2 ry M F(r)[p(r)-p], 



(39) 



where p = $o/B is the mean density. For the Fourier- 
filtered methods 

£cou P © 1 j d 2 r ^^(Q^^expCiQ^rMr) 

S> h x L y J2u(Q»)\PQj 2 - (40) 

For the full method, we have o(r) as a histogram avail- 
able, and we can integrate ( p9[ ) numerically. 

The top plot of Fig. [To] shows how the latent heat varies 
as a function of field. We have shown the jump in inter- 
layer coupling energy normalized by the latent heat in the 
bottom part to demonstrate the contribution of the inter- 
layer coupling to the latent heat. This plot shows that 
the substrate contribution to the latent heat dominates 
at low fields, and becomes less and less important towards 
high fields. 

Fig. [ll] shows the entropy jump across the transition, 
ASp = L p /T m , as a function of field. We find that AS'p 
monotonically decreases with increasing field, as the sys- 
tem approaches the 2D regime. An important issue is the 
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crossover to the 2D melting regime at very large B. Two 
melting scenarios are possible in two dimensions: a usual 
1st order melting and continuous dislocaticp«iiiediated 
melting via the intermediate hexatic phase.EZIu In the 
first case A.S P has to approach a finite value at B — > oo 
and in the second case it should vanish. Early simula- 
tions for a relatively small number (< 500) of logarith- 
mically interacting particles suggested a 1st order phase 
transition.cjL§E3 However, it is known that to resolve 
a continuous melting transition in two dimensions, very 
large systems are required (see, e.g., Ref. [}£]). Therefore, 
the nature of the melting transition of 2D particles with 
logarithmic interactions is an open issue. Resolving this 
issue is beyond the scope of this paper. 

At low fields the entropy weakly diverges for B — ► 0. 
We understand this as follows: the possible configura- 
tions scale as ~ A/t; 2 , where £ 2 is the size of a pan- 
cake, and A is the space it can occupy. For the solid 
state close to the transition the reduced configuration 
space is ~ {u 2 }/£, 2 , because the pancake is confined to an 
area A ~ (u 2 ). In the liquid the reduced configuration 
space grows to ~ a-o/S, 2 , where ao is the average spac- 
ing between pancakes. We get thus an entropy difference 
AS P ~ ln(ag/(u 2 )). Since (u 2 ) approaches at low fields 
a finite field-independent value of the order of A 2 , this 
explains the observed divergence of AS'p at B — ► 0. 

For a precise comparison with the experimentally ex- 
tracted latent heat of vortex-lattice meltingQ one should 
be careful to include the temperature dependence of A, 
which was shown in Ref. ^0| to give extra terms in the 
observable entropy jump. 



cake fluctuation width (it 2 ) and showed how it varies as 
a function of temperature: the variation is significantly 
non-linear below the melting transition, as predicted in 
Ref. |2^. We also calculate the entropy jump across the 
melting transition, which diverges weakly towards small 
fields and large melting temperatures. 

While we have found a satisfying agreement between 
our results here and the earlier approximate work of 
Ref. |2^, the true motivation of this project are the pos- 
sible extensions that can be studied. There is now the 
exciting prospect to study this pancake vortex system in 
the presence of pinning disorder, nThis has been a con- 
troversial topic in recent yearsc2lc3cj'a that our method 
should bring some clarity to. 

Our results cannot be directly compared with experi- 
ments in available layered superconductors because even 
in the most anisotropic BSCCO the Josephson coupling 
is not negligible. However, the position of the melting 
line without the Josephson coupling provides a conve- 
nient reference allowing one to understand the role of the 
Josephson coupling in stabilizing the crystalline phase. 
In particular, it seems that at low fields even a very small 
Josephson coupling such as in BSCCO gives a large up- 
ward shift to the melting temperature. We also note that 
it is possible to suppress the effective Josephson coupling 
by applying a strong in-plane field, as was done in Ref. |45j, 
bringing the melting line closer to the "Josephson-free" 
location. By extending the model to use several layers, it 
is possible to include Josephson coupling between them, 
which would realistically describe an anisotropic layered 
high-temperature superconductor. 



IV. CONCLUSIONS 
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APPENDIX A: CORRELATION CORRECTION 
TO FREE ENERGY 

The correlation correction to the pancake energy, ne- 
glected within the mean-field approach, is given by 

5£ = - Jd 2 r dV £7(r - r', n - n') 5 Pn (r) 5p n > (r') 

(Al) 

The correction to the free energy due to this term is given 
up to second order by 

SFm{5£) -± ^ 

where (...)q implies the mean-field averaging. Substi- 
tuting Eq. (Al) in the last equation and noting that 



{5£ 2 ) ~{5£)l 
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(SS) — we derive for the free energy correction per sity correlation function inside one layer. In the next step 
pancake, bf = aPSF '/ '(L x L y N), 



ATL T L 



J d 2 r d 2 r' J d 2 n d 2 r[ U(r-r',n) 

x U (n - r[,n) (bp (r) bp (n)) (dp (r') bp (r' 1 )) 

where N is the total number of layers, L x L y is the layer 
area, and (bp (r) bp (i"i)) = (<5p„ (r) bp n (ri)) is the den- 

I 



we introduce notation for the sum 

for which, using the mixed representation for 
the interlayer magnetic interaction, U(k±,n) = 

_2^ep n V j.; we obtain the formula 



W(r,ri) 



3 ,-2 



/j2, i2 



exp (ik-r + iki-ri) 



27T 2 Aj " fe 2 fc 2^ + A 2 /c 2 v / 1 + A 2 fc 2 ( + A 2 fc 2 + ^1 + A 2fc2 

allowing us to represent W(r, ri) in a scaling form 

W(r,r 1 ) = ^«;(r/A,r 1 /A). 

By also using a scaling representation for the in-plane density correlation function, (bp (r) bp (i"i)) = \h(r/a, T\/a), 
we derive the scaling representation for the free energy correction 



6f = - S£o ^A G 



a 2s£q 



A' T 



where 



a 


2se " 


1 / 


A ' 


T 


L x L y a J 



l —, [ d 2 rdV d 2 r 1 d 2 r[w 
V 6 J 



r — r' ri — r- 



A ' A 



- I h(r/a, ti/a)h{r' /a,r' x /a) 



(A2) 



(A3) 



r 



is a dimcnsionless function of the order unity. The mean- 
field free energy per pancake, /mf, has the scaling prop- 
erty, 



MF 



SEqG 



a 2seq 



Therefore the free energy correction due to interlayer cor- 
relations ( A2) is smaller than the main term by the factor 
s/A. In particular, the correlation correction shifts the 
melting temperature up as 



bT m = 



bhq -Sf c 
AS 



(A4) 



where bfn q (bf CT ) is the correlation correction to the liq- 
uid (crystal) free energy at the melting point and AS is 
the melting entropy jump. In principle, the mean- field 
simulations allow to compute the correlation correction 
and the corresponding shift of the melting temperature. 
However this computation includes the num erica l evalu- 
ation of a quite cumbersome integral in Eq. (A3). 
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